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ABSTRACT 

In this paper, we analyze the structure of RRab star light curves using Principal Com- 
ponent Analysis. We find this is a very efficient way to describe many aspects of RRab 
light curve structure: in many cases, a Principal Component fit with 9 parameters can 
describe a RRab light curve including bumps whereas a 17 parameter Fourier fit is 
needed. As a consequence we show show statistically why the amplitude is also a good 
summary of the structure of these RR Lyrae light curves. We also use our analysis 
to derive an empirical relation relating absolute magnitude to light curve structure. 
In comparing this formula to those derived from exactly the same dataset but us- 
ing Fourier parameters, we find that the Principal Component Analysis approach has 
distinct advantages. These advantages are, firstly, that the errors on the coefficients 
multiplying the fitted parameters in such formulae are much smaller, and secondly, 
that the correlation between the Principal Components is significantly smaller than 
the correlation between Fourier amplitudes. These two factors lead to reduced formal 
errors, in some cases estimated to be a factor of 2, on the eventual fitted value of the 
absolute magnitude. This technique will prove very useful in the analysis of data from 
existing and new large scale survey projects concerning variable stars. 
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. Kanbur et al (2002), Hendry et al (1999), Tanvir et al (2004) 
' introduced the use of Principal Component Analysis (PCA) 
in studying Cepheid light curves. They showed that a major 
advantage of such an approach over the traditional Fourier 
method is that it is much more efficient: an adequate Fourier 
description requires, at best, a fourth order fit or 9 parame- 
ters, whilst a PCA analysis requires only 3 or 4 parameters 
with as much as 81% of the variation in light curve struc- 
ture being explained by the first parameter. Later, Leonard 
et al (2003) used the PCA approach to create Cepheid light 
curve templates to estimate periods and mean magnitudes 
for HST observed Cepheids. The purpose of this paper is to 
apply the PCA technique to the study of RR Lyrae light 



1 INTRODUCTION 



Walker (2001, hereafter KW). These data consist of 383 
RRab stars with well observed V band light curves in 20 
different globular clusters. KW performed a Fourier fit to 
these data, which, in some cases, is of order 15. Details con- 
cerning the data can be found in KW. The data we work 
with in this paper is this Fourier fit to the magnitudes and 
we assume that the Fourier parameters published by KW 
are an accurate fit to the actual light curves. We start with 
the data in the form used in KW: a list of the mean mag- 
nitude, period and Fourier parameters for the V band light 
curve. The light curve can thus be reconstructed using an 
expression of the form 



curves. 




(1) 



The mathematical formulation and error characteristics 
of PCA are given in K02 and will only be summarized here. 



where Aq is the mean magnitude, u — 27r/P, P the pe- 
riod, Ak,4>k the Fourier parameters given in KW. These light 
curves are then rephased so that maximum light occurs at 
phase and then rewritten as 



The data used in this study were kindly supplied by Ko- 
vacs (2002 private communication) and used in Kovacs and 



2 DATA 




(2) 
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The Uk , bk are the light curve characteristics entering into 
the PCA analysis (K02). We then solve equation (4) of K02, 
either after, or before removing an average term from the 
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Fourier coefficients in equation (2). Witti PC A, the fight 
curve is written as a sum of "efementary" fight curves, 

V{t) = PCAl.Li{t) + PCA2.L2{t) + PCAS.Lsit) + , (3) 

where V(t) is the magnitude at time t, PC Al, PC A2.. etc. 
are the PCA coefficients and the Li{t),i = 1,2,3... are the 
elementary fight curves at phase or time t. These efementary 
fight curves arc not a priori given, t)ut arc estimated from 
the dataset in question. Each star has associated with it a 
set of coefficients PCAl, PCA2, ... and these can be pfotted 
against period just as the Fourier parameters in equation 
(1) are pfotted against period. We afso note that the PCA 
rcsufts arc achieved as a resuft of tfic anafysis of the entire 
dataset of 383 stars whereas the Fourier method produces 
resufts for stars individuaffy. This feature of PCA is partic- 
ufarfy usefuf when performing an ensemble analysis of large 
numbers of stars obtained from projects such as OGLE, MA- 
CHO and GAIA. 



3 RESULTS 

Sofviiig equation (4) of K02 yields tlic Principal Component 
scores and the amount of variation carried by each compo- 
nent. What we mean by this is the following: if we carry 
out an N*'' order PCA fit, then PCA wiU assume that all 
the variation in the dataset is described by N components 
and simply scale the variation carried by each component ac- 
cordingly. Table 1 shows this " amount of variation" quantity 
with and without the average term removed. We see that in 
the case when we do not remove the average term the first 
PC explains as much as 97% of the variation in the light 
curve structure. In the case when we do remove the average 
term from the Fourier coefficients, the first PCA coefficient 
expfains as much 81 percent of the variation in fight curve 
structure. In either case, the first four components expfain 
more than 99.99% of the variation. 

Figures 1 and 2 show some representative fight curves 
from our RRab dataset. In each panel of these two figures, 
the solid line is the Fourier decomposition of order 15 (that 
is 31 parameters) used by KW, whilst the dashed line is a 
PCA generated light curve of order 14 (that is 15 parame- 
ters). Straightforward light curves such as the one given in 
the bottom and top left panels of figures 1 and 2 respectively 
are easily reproduced by our method. The top left panel of 
figure 1 provides an example of an RRab light curve with a 
dip and sharp rise at a phase around 0.8. This is well repro- 
duced by PCA. It could be argued that PCA does not do 
as well as Fourier in mimicking this feature, for example, in 
tlic bottom right panel of figure 2. However, the difference 
in tfic peak magnitudes at a phase of around 0.8 is of the 
order of 0.02mags. It is afso important to remember that the 
PCA method is an ensembfe method and anafyzes aff stars 
in a dataset simuftanoousfy. With Fourier, it is possibfe to 
taifor a decomposition to one particufar star. This differ- 
ence can be seen cither as a positive or negative point about 
either technique. Given this, we contend that PCA does re- 
markably well in describing the full light curve morphology 
of RRab stars. On the other hand, the Fourier curve in the 
bottom left panel of figure 2 at this phase is not £is smooth 
as the PCA curve. 

In fact the PCA curves do not change much after about 



8 PCA parameters. Even though table 1 implies that the 
higher order PCA eigenvalues are small, we feel justified in 
carrying out such a high order PCA fit because its only 
after about 8 PCA components that the fitted light curve 
assumes a stable shape. The left panel of figure 3 displays 
an eighth order PCA fit (9 parameters, dashed line) and 
a fourth order Fourier fit (9 parameters, solid line). The 
Fourier curve still has some numerical wiggles whilst the 
PCA curve is smoother. In addition, the two curves disagree 
at maximum light. The right panel of figure 3 shows, for the 
same star, the same order PCA curve as the left panel and an 
eighth order Fourier fit (17 parameters). Now the two light 
curves agree very well. Note that in portraying the PCA 
and Fourier fits of reduced order in this figure, we simply 
truncated the original representations to the required level. 

We suggest that figures 1-3 and table 1 provide strong 
evidence that PCA is an efficient way to describe RRab light 
curve structure without compromising on what light curve 
features are captured by this description. 

Figures 4-6 display plots of the first three PC scores 
plotted against log period for our sample. The errors asso- 
ciated with these PCA scores are discussed in section 4 of 
K02 and given in equation 6 of that section. The orthogonal 
nature of tlicsc scores may well provide insight into the phys- 
ical processes causing observable features in the light curve 
structure. A detailed study of these plots, in conjunction 
with theoretical models, is left for a future paper. 

Figure 7 graphs V band amplitude against the first PCA 
coefficient (after averaging). We see a very tight correfation. 
Since tabfe 1 implies that PCAl explains about 81% of the 
variation in light curve structure, figure 6 shows that the 
amplitude is a good descriptor of RRab light curve shape, 
at least for the data considered in this paper. Although the 
Fourier amplitudes are also correlated with amplitude, with 
PCA, we can quantify, very easily, the amount of variation 
described by each PCA component. This has implications 
for both modeling and observation. On the modeling side, a 
computer code that can reproduce the observed amplitude 
at the correct period, will also do a good job of reproduc- 
ing the light curve structure. On the observational side, this 
provides insight into why we can use the amplitude, rather 
than a full blown PCA or Fourier analysis, to study the gen- 
eral trends of light curve structure. This is why comparing 
theoretical and observational RRab light curves on period- 
amplitude diagrams works reasonably well, though we cau- 
tion that a careful analysis should consider the finer details 
of light curve structure. 

Figures 6 and 7 display plots of the first two PCA co- 
efficients and Fourier amplitudes, respectively, for our data, 
plotted against each other. Whilst Ai and A2 are correlated 
with each other, PCAl and PCA2 are not, by construc- 
tion. A similar situation would occur had we plotted Ai or 
A2 against A3. This is another advantage of PCA analysis 
of variable star light curves: the different PCA components 
are orthogonal to each other. A practical advantage of this 
feature is outlined in the next section. 

4 LIGHT CURVE LUMINOSITY RELATIONS 

A major goal of stellar pulsation studies is to find formu- 
lae linking global stellar parameters such as luminosity or 
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phase phase 



Figure 1. Light curve reproduction using Fourier (solid lines) and PC A (dashed lines) methods 



metallicity to structural light curve properties. If we are in- 
terested in the V band magnitude, then we can write, 

Mv = f{lightcurvestructure), 

where, since we do not know the function /, we try to esti- 
mate it empirically. Two different approaches to quantifying 
light curve structure will, in general, yield different formu- 
lations of the function /, but if there does exist a true un- 
derlying function /, then both methods should give similar 
answers for M^, given the same input data. With a Fourier 
based method, the function / is related to the Fourier am- 
plitudes and phases, Ak,<j)ki, usually with a linear relation. 
With a PCA approach, we use the PCA scores plotted in fig- 
ures 2-4. Hence a PCA relation, though also linear, will be 
different. The nature of PCA implies that the error structure 
in such formulae will bo simpler and wc quantify this below. 
Both formulations should, of course, give similar numbers 
for the final estimated value of the physical parameter in 
question, in this case. My. 



KW used the Fourier method and found relations of the 
form. 

My = const. - 1.82 log P - 0.805Ai , (4) 
and. 

My = const. - 1.876 logP - 1.158Ai -|- 0.821^3- (5) 

We note that these relations were obtained through an 
iterative procedure whereby outliers were removed and the 
relations re-fitted (Kovacs 2004). In this paper, we use the 
PCA method, but also, we use the entire datasct C men- 
tioned in KW, consisting of 383 stars, and fit the relations 
just once. We do not remove any outliers. This may be why 
wc obtain slightly different versions of the fit using Fourier 
parameters than that published in KW. For ease of com- 
parison, wc include in tabic 2 results obtained using both 
PCA and Fourier parameters. This table gives the name for 
the relation, the independent variables considered and co- 
efficients together with their standard errors. The value of 
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Figure 2. Light curve reproduction using Fourier (solid lines) and PCA (dashed lines) methods 



chi-squared in the table is defined as 

k=N 

J2iMv-M,f/iN-p), (6) 

fc=i 

where Mv is the fitted value of Mv and A^, p are the number 
of stars and parameters respectively in the fit. An examina^ 
tion of this table strongly suggests that 

• 1) Similar relations to equations (4) and (5) between 
M„ and the PCA coefficients exist. 

• 2) Wc can use an F test (Weisberg 1980) to test for 
the significance of adding a second and then a third PCA 
parameter to the regression. The F statistic we use is 

{RSSnh — RSSAH)/{dfNH — dfAn) 

RSSAn/dfAH, 

where RSSnh, RSSah are the residual sum of squares un- 
der the null and alternate (NH and AH) hypothesis respec- 
tively. Similarly, d/jv/f and dfAH are the degrees of freedom 



under these two hypotheses. For this problem, the null hy- 
pothesis is that the model with the smaller number of pa- 
rameters is sufficient whilst the alternative hypothesis is that 
the model with the greater number of parameters is required. 

Under the assumption of normality of errors, equation (7) 
is distributed as an ^(d/^^jj— df^jij^),^/^^^^ ? 

(Weisberg, 1980, p. 

88). Applying this F test implies firstly, that adding the 
first parameter PCAl is a significant addition to logP and 
secondly, that adding a second and third parameter, PCA2 
and PCA3 arc also highly significant with a p value loss 
than 0.0004. In the case of Fourier parameters, adding the 
Ai parameter to log P is highly significant and adding the 
As parameter to this is also highly significant. However, a 
formula involving (log P, Ai,A2) has a p value of 0.0058 and 
a formula involving all 3 Fourier amplitudes and log P is not 
a significant addition to a formula involving (log P, Al, A3). 

• 3) The standard deviation of the fits given in the last 
column is generally slightly higher for the PCA case, when 
considering similar numbers of parameters. This is perhaps 



Principal Component Analysis of RR Lyrae light curves 5 





Figure 3. Light curve reproduction using Fourier (solid lines) and PCA (daslied lines) methods. The left panel is a fourth order (9 
parameters) Fourier fit and an eight order PCA (9 parameters) fit. The right panel is an eight order (17 parameters) Fourier fit and an 
eight order PCA (9 parameters) fit. ( 
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Figure 4. Plot of first Principal Component against log period. 



Figure 5. Plot of second Principal Component against log period. 



caused by the fact that the different PCA components carry 
orthogonal sets of information. 

• 4) The errors on the coefficients in the PCA fits are 
always significantly smaller. This is an important point when 
we evaluate the errors on the final fitted value of the absolute 
magnitude. 

• 5) If wc write the absolute magnitude as a function of 

parameters, xi,X2, --jXn, 



then the error on the absolute magnitude is given by, 
fT^(Mv + const.) = V a^{xk){^) + 

^-^ OXk 

k=l 

t 



(9) 



My + const. = f{xi,X2, xn), 



(8) 



As table 2 indicates, a^{xk) is always smaller when the Xk 
arc PCA coefficients rather than Fourier amplitudes. Figure 
8 and 9 portray graphs of PCAl vs PCA2 and Ai verses A2 
respectively. We note that pija{xi)a{xj) = (T^{xi,Xj). Ta- 
ble 3 presents sample correlation and covariance coefficients 
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Figure 8. Plot of first Principal Component against second Prin- 
cipal Component. 



V amplitude 



Figure 7. Plot of V band amplitude against the first PCA co- 
efficient. 



between the period and PCA parameters and period and 
Fourier parameters. Table 3, and figures 6 and 7 demonstrate 
that the correlation coefficient amongst any pair of PCA co- 
efficients is smaller than between any pair of Fourier coeffi- 
cients. Hence the error on the fitted value of M„, (j^(M„), has 
to be smaller when using a PCA based formula. We can use 
table 3 and equation (9) to formally calculate the error on 
M„ -|- const. Table 4 presents these results. The label in the 
top row of this table (PI, Fl, etc.,) refers to the appropriate 
relation in table 2. Wc sec clearly that the PCA formulae do 
better than their Fourier counterparts with a similar number 
of parameters. When we consider the (log P, PCI, PC2) and 
(log P, Ai, As) variables, then the "error advantage" using a 



Figure 9. Plot of first Fourier amplitude against second Fourier 
amplitude. 



PCA based method is a factor of two. This occurs not just 
because the PCA coefficients arc orthogonal to each other, 
but also because the errors on the coefficients in a PCA 
based formula are significantly smaller than in the Fourier 

case. 

Figure 10 displays a plot of the predicted ab- 
solute magnitudes obtained using a two parameter 
(logP, ^1,^3) Fourier fit and the three parameter 
l\ogP,PCAl,PCA2,PCA2,) PCA fit. The two approaches 
are displaced from each other because we do not consider the 
constants in this study. Disregarding this, it can be seen that 
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Figure 10. Plot of fitted My + const values when using Fourier 
and PCA methods. 
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Figure 11. Plot of absolute magnitude difference versus frac- 
tional change in light curve parameters for Fourier (open squares) 
and PCA (closed squares). 



the slope of this plot is 1: hence the two methods produce 
similar relative absolute magnitudes. 



5 CONCLUSION 

We have shown that the method of PCA can be used to 
study RR Lyrae light curves. It has distinct advantages over 
a Fourier approach because 



• a) It is a more efficient way to characterize structure 
since fewer parameters are needed. A typical Fourier fit re- 
quires 17 parameters whereas a PCA fit may only need 9. 

• b) Using the PCA approach, we see clearly why the 
amplitude is a good descriptor of RRab light curve shape. 

• c) The different PCA components are orthogonal to 
each other whereas the Fourier amplitudes are highly corre- 
lated with each other. This leads to relations linking light 
curve structure to absolute magnitude using PCA having 
coefficients with smaller errors and leading to more accu- 
rate estimates of absolute magnitudes. This can reduce the 
formal error, in some cases, by a factor of 2. 

In the present formulation of our PCA approach, the 
input data is a Fourier analysis. If these input data, that 
is the Fourier decompositions, contain significant observar- 
tional errors, the error bars on the resulting Principal Com- 
ponents will be larger. Neither the PCA or Fourier approach 
can compensate fully for noisy data. In this sense, the sen- 
sitivity of PCA to noisy data should be similar to Fourier, 
though the fact that PCA is an ensemble approach in which 
we initially remove an average term does guard against in- 
dividual points having too much undue influence. As an ex- 
ample, table 4 of KW gives 17 outliers (in terms of their 
Fourier parameters), which KW removed in their analysis 
relating absolute magnitude to Fourier parameters. We do 
not remove these outliers, yet, in terms of the final fitted 
magnitudes presented in figure 10, PCA and Fourier pro- 
duce very similar results. Further, even with the inclusion of 
these 17 stars, the PCA method still produces PCA coeffi- 
cients with smaller errors as given in tables 2 and 3. Kanbur 
et al (2002) discuss in detail the nice error properties of the 
PCA method as applied to variable stars and give a recipe 
with which to calculate errors on PCA coefiicients. Their 
figure 2, albeit for Cepheids, displays error bars on these 
coefiicients. We see that even with noisy data, the progres- 
sion of PCA parameters with period is preserved, though of 
course, the error bars on the PCA coefficients are larger. 

Ngeow et al (2003) developed a simulated annealing 
method which can reduce numerical wiggles in Fourier de- 
composition of sparse data. Ngeow et al (2003) give specific 
examples of how such an approach improves Fourier techin- 
ues using OGLE LMC Cepheids. A similar result will hold 
true for RR Lyraes. Hence this annealing technique cou- 
ple with a Principal Component analysis should prove very 
useful when dealing with noisy RR Lyrae data and will be 
treated in detail in a subsequent paper. 

Our PCA results are based on a sample of 383 stars in 
globular clusters. How transferable are our results and how 
can our results be used to obtain PC coefficients for a new 
RR Lyrae light curve which appears to be normal (ie no 
signs of Blazhko effects etc.)? 

Our results are transferable to the extent that the orig- 
inal 383 stars are a good representation of the entire popu- 
lation of RRab stars, including variation in metallicity and 
differences between field and cluster variables. Given this 
caveat, we suggest two methods to reproduce the light curve 
of a new RRab star. Firstly, it is straightforward to include 
the new star in the PCA analysis with the existing dataset. 
This is our recommended approach and preserves the "en- 
semble analysis" property of our PCA method. Our second 
method will be the subject of future paper but briefly it is 
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this. We fit the progression of the PCA coofEcicnts with pe- 
riod, such as given in figures 4 and 5, with simple polynomial 
functions. As an aside, we remaxk that figure 4 contains sig- 
nificant scatter, perhaps associated with metallicity, so that 
it would bo best to include metallicity in such polynomial 
fits. For a new star, we then guess its period and read oflF, for 
that period, the value of the PCA coefficients. Equation (3) 
then allows us to generate the light curve. We iterate this 
until a specified error criterion is satisfied. We can then use 
existing formulae relating absoltc magnitude to light curve 
structure as defined by PCA. This PCA template approach 
has been used, with considerable success, in analysing HST 
Cepheid data (Leonard ot al 2003). 

We note from table 2 that the chi square on the fitted 
relations arc similar for PCA and Fourier. Does this mean 
that despite the smaller formal errors with PCA, both meth- 
ods' ability to predict RRab absolte magnitudes is limited by 
the intrinsic properties of RRab stars themselves? To some 
extent this is true. Jurcsik et al (2004), in analysing accurate 
data for 100 RRab stars in M3, show that for some 16 stars, 
amongst which there exist some pairs whose absolute mean 
magnitudes difi^er by about 0.05 mags (the accuracy of the 
photometry is about 0.02mags), the Fourier parameters and 
periods are very similar. That is, an empirical method re- 
lating absolute magnitude to period and Fourier parameters 
in one waveband could not distinguish between these stars. 
Since, as Jurcsik et al (2004) point out, their data contains 
a small range of both mass and metallicity, temperature is 
the only other variable, it may be the case that multiwave- 
length information is needed. It is worthwhile to investigate 
how PCA fares with this dataset. Here we give an outline 
that suggests that PCA can be more efficient at extracting 
information from the light curve. 

For the sixteen stars which had differing absolute mag- 
nitudes but very similar Fourier parameters, we can perform 
the following procedure: for every pair, j ^ k, we calculate 

(ol(i) - ol(fe))/ol(fc) -I- (o2(j) - a2{k))/a2{k) 
+{a3ij) - o3(fe))/a3(fc) = diffl, 
(pcal(j) — pcal{k)) /pcal{k)+ 

{pcal{j)-pca2{k))/pca2ik) = diff2, 

and 

{vmean{j) — vmean{k)) = diff3, 

where al{j),a2(j),a3{j) are the Fourier amplitudes and 
pcal{j) , pca2{j) arc the PCA coefficients and vmean(j) arc 
the mean magnitudes. In the above, we always take the ab- 
solute value of the differences. We need to take fractional 
changes because the Fourier amplitudes and PCA coeffi- 
cients have different ranges. We now plot diff3 against diffl 
and diff2. This is presented in figure 11, where the open 
squares are diffl and the closed squares are diff2. We see that 
with PCA (closed squares), the differences between light 
curve structure parameters are greater than with Fourier 
(open squares). This could imply that PCA can be more 
efficient though the limitations associated with using a sin- 
gle waveband are still present. A more rigorous, quantitative 
discussion of this, in a Fisher information sense, will be given 
in a future paper. 



Table 1. Percentage of variation explained by PC components 



PCI PC2 PCS PC4 PCS PC6 

without average 81.4 7.8 5.7 2.3 0.74 0.57 
with average 96.9 1.9 0.55 0.25 0.07 0.006 



In other future work we plan to investigate the appli- 
cability of this method to light curve structure-metallicity 
relations, RRc stars and a comparison of observed and the- 
oretical light curves using PCA. 
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Table 2. Light curve luminosity relation using PCA and Fourier methods. 





logP 


first 


second 


third 


chisquare 


PCA 


PO 


-1.134±0.059 








0.00321 


PI 


-1.550±0.082 


0.269±0.038 






0.00283 


P2 


-1.609±0.082 


0.290±0.038 


0.291±0.082 




0.00274 


P3 


-1.744±0.088 


0.329±0.039 




-0.539±0.107 


0.0027 


P4 


-1.829±0.088 


0.359±0.039 


0.336±0.079 


-0.583±0.105 


0.00253 


Fourier 


Fl 


-1.677±0.083 


-0.472±0.054 






0.00266 


F2 


-1.700±0.082 


-0.726±0.092 




0.613±0.179 


0.00258 


F3 


-1.740±0.085 


-0.758±0.116 


0.536±0.193 




0.00261 


F4 


-1.720±0.085 


-0.790±0.117 


0.215±0.243 


0.490±0.227 


0.00258 



Table 3. Sample correlation and covariance coefficients between period, PCA and Fourier coefiicients 





log P, PC Al 


log P, PCA2 


log P, PC A3 


PCA1,PCA2 


PCA2, PC A3 


PCAl, PC A3 


correlation 


0.631 


0.099 


-0.299 


< 10-8 


< 10-6 


< 10-^ 


covajriance 


0.0038 


0.0002 


-0.0006 


< io-« 


< io-« 


< io-« 




logP.Ai 


logP.Aa 


log P, A3 


Ai,A2 


M,A3 


Ai,A3 


correlation 


-0.655 


-0.529 


-0.562 


0.926 


0.931 


0.902 


covariance 


-0.0028 


-0.0012 


-0.0011 


0.0030 


0.0013 


0.0024 



Table 4. Formal error on My + const, for PCA and Fourier relations 



PI P2 P3 P4 Fl F2 F3 F4 

0.0139 0.0142 0.0216 0.0240 0.0156 0.0313 0.0394 0.0311 



